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Quantum Monte-Carlo (QMC) simulations involving fermions have the notorious sign problem. 
Some well-known exceptions of the auxiliary field QMC algorithm rely on the factorizibility of the 
fermion determinant. Recently, a fermionic QMC algorithm 1J has been found in which the fermion 
determinant may not necessarily factorizable, but can instead be expressed as a product of complex 
conjugate pairs of eigenvalues, thus eliminating the sign problem for a much wider class of models. 
In this paper, we present general conditions for the applicability of this algorithm and point out that 
it is deeply related to the time reversal symmetry of the fermion matrix. We apply this method to 
various models of strongly correlated systems at all doping levels and lattice geometries, and show 
that many novel phases can be simulated without the sign problem. 

PACS numbers: 02.70.Ss, 71.10.Fd 
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I. INTRODUCTION 



Understanding the physics of strongly correlated many 
body systems is a main focus of condensed matter physics 
today. However, most models with strong interactions 
can not be solved exactly except in one dimension. 
Presently, there are no systematic non-perturbative ana- 
lytic methods which work in higher dimensions. Largely 
because of this reason, numerical simulations such as ex- 
act diagonalization (ED), density-matrix renormalization 
group (DMRG), quantum Monte-Carlo (QMC) are ex- 
tensively performed to study strongly correlated systems. 
However, each of the numerical methods has its own lim- 
itations. The ED can only be performed in a very small 
sample size, and the DMRG method is largely restricted 
to one-dimensional systems. In contrast, QMC simula- 
tion is the only systematic and scalable method with suf- 
ficient numerical accuracy for higher dimensional prob- 
lems. However, QMC also has the notorious fermion sign 
problem which makes low temperature properties inac- 
cessible. 

In lattice systems, a particular version of QMC uses 
the auxiliary field method introduced by Blankenbecler, 
Scalapino and Sugar |2(, with fruitful results. Because 
one can not directly sample the fermionic Grassmann 
fields, the standard process is to perform a Hubbard- 
Stratonovich (HS) transformation to decouple the four- 
fermion interaction terms, and then to integrate out the 
fermions 0. The resulting fermion functional determi- 
nant works as the statistical weight for sampling the aux- 
iliary fields. However, generally speaking the fermion de- 
terminant may not be positive, and can even be complex 
in some cases. The sign or the phase of the fermion deter- 
minants can lead to dramatic cancellations which makes 
statistical errors to scale exponentially as the inverse of 
the temperature and size of the system. This notorious 
sign problem is the major obstacle in applying QMC to 
fermionic systems. A successful solution to the sign prob- 
lem would obviously lead to great advances in quantum 
many body physics. 



There are a few exceptions where the sign problem 
is absent, such as the negative U Hubbard model and 
the positive U Hubbard model in a bipartite lattice at 
the half-filling In both cases, the fermion determi- 
nant after the HS decomposition can be factorized into 
to two real parts with the same sign. It is therefore 
positive-definite. Unfortunately, general fermion deter- 
minants may not be factorizable for more complicated 
models and the majority of models do have the sign prob- 
lem. In recent years, several other algorithms have been 
roposed which partially solves the minus sign problem 



Recently, it has been shown that the minus sign prob- 
lem can be eliminated without relying on the factorizibil- 
ity of the fermion determinant, therefore, a broader class 
of models can be simulated by the QMC algorithm p|. 
The fermion determinant can always be expressed as a 
product of its eigenvalues; under certain conditions, the 
eigenvalues of the fermion determinant always appear in 
complex conjugate pairs, thus making the fermion deter- 
minant positive definite. In this article, we shall show 
that the property of conjugate eigenvalue pairs follows 
from the time reversal symmetry of the HS decoupled 
Hamiltonian, and can be viewed as a generalization of the 
Kramers theorem in quantum mechanics. We shall call 
this method the T -invariant decomposition (time rever- 
sal invariant decomposition) . This method does not lead 
to any improvement for the single band Hubbard model, 
but significantly extends the applicability of the QMC 
to multi-band, multi-layer or higher spin models. This 
algorithm is particularly useful for Hubbard models with 
higher spins, which can be accurately realized in systems 
of cold atoms. Recently, Assaad et al. 9] applied the 
QMC to generalized Hubbard models with more bands. 
Imposing the factorizibility condition of the fermion de- 
terminant, they found that they could extend the pa- 
rameter regime for QMC free of the sign problem only 
by scarifying the spin rotational invariance. However, 
applying our method of T-invariant decomposition with- 
out requiring factorizibility, we shall show that multi- 
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band or higher spin Hubbard models can be simulated 
for an extended parameter regime without scarifying the 
spin rotational invariance. This QMC algorithm based 
on T-invariant decomposition has been recently applied 
to conclusively demonstrate the stag gere d current carry- 
ing ground state in a bi-layer model |lfl |. 

The rest of this article is outlined as follows: In Sec. 
ITT1 the sign problem for the spin 1/2 Hubbard model is 
reviewed. In Sec. IIIII we prove the fundamental theo- 
rem of T-invariant decomposition and show the absence 
of the sign problem. In Sec. IIVI we employ the algo- 
rithm to the spin 3/2 Hubbard model and the general- 
ized arbitrary spin n — 1/2 fermionic Hubbard model. In 
Sec. E] we apply it to a bi-layer model introduced by 
Scalapino, Zhang and Hanke[ll|, which can be mapped 
into the spin 3/2 Hubbard model. In Sec. IVII we discuss 
the algorithm in the model Hamiltonians with bond in- 
teractions and various exotic phases. Final conclusions 
are presented in Sec. IVIII 



II. THE SIGN PROBLEM IN THE SPIN 1/2 
HUBBARD MODEL 



In this section, we review the sign problem in the spin 
1/2 Hubbard model and interpret its absence in the neg- 
ative U case as due to its time reversal properties of the 



HS decomposition. The Hubbard model on the lattice is 
commonly defined as 

H = -t^2(4 a c ja + h.c.) - flinty 

+ uJ2int(i)-l)(ni{i)-l), (1) 

i 

with t the hopping integral, \i the chemical potential, 
a =T,j> n a (i) = v ia c ia and n(i) = n T (i) + nj.(z). At 
half-filling and on a bipartite lattice, the particle-hole 
symmetry ensures that fi = 0. 

To perform the QMC simulation, we first need to 
decouple the 4-fermion interaction terms using the HS 
transformations by the Gaussian integral: 

1 f 1 

exp(-A 2 ) = V2tt J cfoexp(--x 2 - xA), 

1 f 1 

exp(--yl 2 ) = V27r / dxexp(--a; 2 - ixA). (2) 

Various HS decoupling schemes are discussed in Ref . ^2 • 
For U < 0, it is convenient to decouple Eq. (Q) in the 
density channel and then integrate out the fermions. The 
resulting partition function is given by 



Z = j DJDc exp{- J\t {c\^-c a + H)} 

= yDn^Dc cxp | --^^ dr^(n(i,r)-l) 2 |cxp|-^ dr (H K + #j(t))} 

= j Dnexpj--^ f dr^O(i,T) - l) 2 } dct{/ + B}, 



(3) 



where n(i,r) is a real HS bose density field. The Note that the matrix kernels hftj aa , and hjj aa , entering 
imaginary-time independent kinetic energy term H K in Eq. © , as well as the I + B matrix itself, are 2 N x 2 N 



and the imaginary-time dependent decoupled interaction matrices, if the lattice system under simulation has N = 

L x x L y sites. In the subsequent discussions, we shall 
simply use the second quantized operators Hk and Hj 
interchangeably with the first quantized matrix kernels 
hx and hi to save some writing, whenever their meanings 
are obvious from the context. 



term Hj(t) can be expressed as 

Hr^- S Tc i h K ,c- • Hr-Vc* h 1 ,c- * 

ij i 



h\ jaa , = Un{i,T)5ij5 aal . 

Here aa , and atT , are defined for both spin compo- 
nents on each site. After integrating out the fermions, 
we obtain 



(4) 



In practice, I + B needs to discretized as 



J_|_ B = I + e -ArHK e -ATH"i(n) e -ATir Ke -ArH<(T i _ 1 )_ 



e -ArH K g-ArfliCn) 



(0) 



where At — (3/1 is the discretized time slice. 
I + B = I + Tvcp{- J dr (h K + hj(T))}. (5) Similarly, at [/ > 0, Eq. JTJ can be decomposed in the 
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spin density channel as 

Z= J DS z exp{-2U J dr^5^(i,r)}det{7 + B}, 

(7) 

with the same expression for B as in Eq. JSJ, but with 
Hi replaced by 

Hj(r) = -2Uj2{ct(rX c^(T)}sS,r). (8) 

i 

It is well known that the spin 1/2 Hubbard model is 
free of the sign problem either for U < or for U > 
at half- filling and in a bipartite lattice 0, The 
usual proof is based on the factorization of the fermion 
determinant as 

det{I + B} =det{7 + S T }det{7 + BJ. (9) 

In the negative U case, the HS decomposition in Eq. @ 
enables such a factorization, and Bj is identical to B± for 
any HS field configurations. Therefore det{7 + B} is the 
square of a real number and thus positive definite. Gen- 
erally speaking, in the positive U case, the HS decomposi- 
tion in Eq. © still enables factorization, but det{/+i?|} 
is different from det{7 + B{\, thus the sign problem ap- 
pears. However, at half-filling and on a bipartite lattice, 
it is possible to change the sign of U while keeping the 
kinetic energy part invariant by a partial particle-hole 
transformation only on spin down particles 

c^->c^, ai -> (-) l c- x , (10) 

then the above algorithm is also applicable. Nevertheless, 
this transformation can not be applied to lattices which 
are not bipartite, or away from the half- filling (/i =/= 0), 
thus the sign problem remains in general. 

Recently, an anisotropic two band model explicitly 
breaking the spin rotational symmetry is also shown to be 
free of the sign problem |9j . The Hamiltonian is defined 
by 

H = -t^2{4 a c ja + h.c.) -fj,^2n a (i) 

ij,0 i,cr 

-- \U\Y,{nx{i)-n 2 {i) + m{i)-n A {i)f, (11) 

i 

where n a {i) = c\(i)c a (i) are particle densities for each 
spin components a = 1,2,3,4. The interaction part can 
be decoupled as 

Z = J L>Scxp{-|[/| J P dTj2s 2 {i,T)} 

rf) 

x exp{- / dT(H + Hj(t))} 
Jo 

H i{t) = ^2(4^1 - 4,2°i,2 + 4,3 C h3 ~ 4a c ^a) S (^ t ) 



This HS decomposition enables the factorization of the 
fermion determinant as 

det{7 + B} = det{I + B} 12 det{I + B} 34 , (12) 

where det{/+i?}i2 and det{/+B}34 for spin components 
1,2 and 3,4 respectively are identical and real. There- 
fore, the fermion determinant is positive in this case as 
well. However, a disadvantage of this model is the ex- 
plicit breaking of the spin rotational symmetry. 

III. FUNDAMENTAL THEOREM OF 
T INVARIANT DECOMPOSITION 

We now show that the condition of factorizibility of 
the fermion determinant is unnecessarily restrictive, and 
a more general condition can be precisely stated. The 
fermion determinant is a product of all the eigenvalues. 
Since I + B involves a time ordered product, it may not 
be Hermitian and the eigenvalues may be complex in gen- 
eral. Because the ensemble of HS field configurations is 
arbitrary, one would naively not expect any special rela- 
tions among the eigenvalues. Surprisingly, the time rever- 
sal symmetry provides an important relationship among 
the eigenvalues. To formulate the fundamental theorem, 
we consider Hk and Hi in the I + B matrix of Eq. ijfjj to 
be the HS decomposed single particle Hamiltonian ma- 
trix derived from a general Hamiltonian, not necessarily 
the s = 1/2 Hubbard model. 

Theorem: If there exists an anti-unitary operator T, 
such that 

THkT- 1 = H K , THff- 1 = H I: T 2 = -1, (13) 

then the eigenvalues of the I + B matrix always appear 
in complex conjugate pairs, i.e., if is an eigenvalue, 
then A* is also an eigenvalue. If Ai is real, it is two- 
fold degenerate. In this case, the fermion determinant is 
positive definite, 

dct(/ + B) =JJ|A j; | 2 > 0. (14) 

i 

Proof: From the condition of the theorem stated in Eq. 
(H3J), it obviously follows that T(I + B)T _1 = (I + B). 
For simplicity, we first consider the case where J + B is 
an n x n dimensional diagonalizable matrix, i.e., there 
exists a non-singular matrix P satisfying 

P-\I + B)P - diag{Ai, A 2 , A„}. (15) 

The n columns of P can be viewed as a set of linearly 
independent state-vectors 

P = {!*!>, |¥ 3 >,.. .,|* n », (16) 

Suppose that is an eigenvector with eigenvalue Aj, 
i.e. (I + B)\^>i) = Ai|^i). Using the anti- unitary prop- 
erty of T, we see that 

(I + B)T\*i) = T(I + B)T- 1 T|1' i ) = X*T\^i). (17) 
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Therefore, T\*f?i) is also an eigenvector, with eigenvalue 
A*. Since T 2 = -1, and \%) are orthogonal to 

each other. This shows that Xi and A* are two different 
eigenvalues, thus the eigenvalues of I + B appear in com- 
plex conjugate pairs as stated in the theorem. If I + B is 
Hermitian, our theorem reduces to Kramer's theorem on 
the time reversal symmetry in quantum mechanics, stat- 
ing that the eigenvalues of / + B are real, and two-fold 
degenerate. 

In the general case, I + B may not be diagonalizable, 
instead it can always be transformed into the Jordan nor- 
mal form as diagonal blocks 



P-\I + B)P = diag{ J l7 J 2 , J fe }, 



(18) 



where P is an n X n non-singular matrix as before and 
Ji is an li x li bi-diagonal matrix as 



A,; 1 



Ji 



V 



1 



\ 



J 



(19) 



The determinant of / 
eigenvalues 



B is still the product of all the 



det(/ + B) = jj(Ai)' 



(20) 



i=l 



As in Eq. 1161 P can be viewed as n linearly independent 
column state- vectors as 



P — {Pi, Pii Pk}> 



(21) 



where each Pi is an n x li matrix containing li column 
state- vectors 



Pi = {|* m+1 ),...,|* m+ii >}, 

For each Jordan block Ji, it satisfies 
(I + B)P =PJi, 



i-1 



(22) 



(23) 



thus among the state- vectors in Pi, \^ m+ i) is the only 
eigenvector with eigenvalule Xi. It is straightforward to 
show that 

(I + B)(TPi) = (TPjj;, (24) 

where (TPi) is defined as 

(TP i ) = {T\* m+1 ), T\* m+h )}, (25) 

and T|\I> m+ i) is the only eigenvector with eigenvalue X* 
in (TPi). Again since l^m+i) and T|\I/ m+ i) are orthog- 
onal to each other, (TPi) contains different state vectors 
from what Pi does. As a result, Ji and J* are different 
Jordan blocks. As before, the Jordan blocks appear in 
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FIG. 1: Distribution of eigenvalues in the complex plane, (a) 
Eigenvalues of a fermion matrix satisfying the conditions of 
our theorem are always paired, (b) Complex eigenvalues of 
a generic real matrix are paired, but real eigenvalues are not 
two-fold degenerate in general, leading to negative determi- 
nants. 



complex conjugate pairs and so do the eigenvalues. This 
completes the proof for the general case of / + B. 

Since the anti-unitary operator T used in our theo- 
rem shares similar properties as the time reversal trans- 
formation in quantum mechanics, we call our method 
T-invariant decomposition. However, it is important to 
emphasize that any anti-unitary operator with the stated 
mathematical properties could work here. In some exam- 
ples we shall discuss, T does not have the explicit physical 
meaning of the time reversal transformation. 

It is also important to point out that the T 2 = — 1 
condition is essential for our theorem. In the case when 
the fermion matrix is real, one can define a trivial anti- 
unitary operator T — C, where C denotes the complex 
conjugation. In this case, if the eigenvalue A^ is com- 
plex, i.e., Xi y^z X*, then A* must also be an eigenvalue. 
However, when is real, it is in general not two- fold de- 
generate, since |^) and T\^>) may not be orthogonal for 
the case of T 2 = 1. In this case, an odd number of neg- 
ative eigenvalues would lead to a negative determinant. 
The distribution of eigenvalues in the complex plane for 
a fermion matrix satisfying the condition of our theorem 
and the eigenvalues of a generic real fermion matrix is 
illustrated in Fig. ^ When the conditions of our theo- 
rem is violated, either the complex conjugate eigenvalue 
pairs collide on the real axis and move off from each other 
along the real axis, or the two-fold degenerate eigenvalues 
move off directly from each other along the real axis. 

A restricted version of our theorem was originally dis- 
cussed in the context of nuclear physics|5j- However, 
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these authors overlooked the case that I + B may not 
be diagonalizable, thus their proof was not complete. In 
addition, our T-transformation is not restricted to the 
physical time-reversal transformation as in Ref.Q, thus 
the theorem applies to a much wider class of models. 

We now illustrate this general theorem for the case of 
the s = 1/2 Hubbard model. For the spin 1/2 system on 
each site, the time-reversal transformation T is defined 
as T(i) = R(i)C, satisfying T 2 (i) = -1 , where 

R = -ia v = ( J - 1 V (26) 

For the entire system, the time-reversal operator is de- 
fined as the direct product T = (Hi (g) R(i))C. The four 
independent fermion bilinears in the particle-hole channel 
can be classified as the particle number n(i) — ^| a '0j,a 

and spin S(i) = ip\ a (u /2) Q , | a^i ll g, which are even and odd 
under the T transformation respectively: 

Tn^T- 1 = n(i), TS^T- 1 = -S(i), (27) 

Now we can understand the absence of sign prob- 
lem in the negative U case as follows. The density 
channel decomposition is T-invariant, namely, T(Hk + 
i?/(r))r" 1 = H K + Hj{t). The conditions of our theo- 
rem is satisfied and the fermion determinant is thus pos- 
itive. For U > 0, the Hamiltonian can be decoupled in 
the density channel at the cost of involving the imaginary 
number i, or decoupled in the spin channel with only real 
numbers. In either cases, while Hk is still even under T, 
Hj is odd. The conditions of our theorem does not apply, 
and the sign problem appears in general. 

For a general interacting fermion model, we can always 
express T = RC, where RR* = — 1 and R* is the complex 
conjugate of R. In many cases, R is purely real, and 
it reduces to R 2 = — 1. The general condition for our 
theorem then reads 

R(H K + H^R- 1 = {H K + H T )*, (28) 

with the unitary matrix R satisfying RR* = — 1 for any 
configurations of the HS field. Again we emphasis that 
the precise form for R in Eq. (|26|l is not necessary. 

While our new method does not lead to any improve- 
ment of the sign problem for the s — 1/2 Hubbard model, 
we shall show now that it significantly improves the QMC 
algorithm for multi-band, multi-layer and higher spin 
models, since the conditions for our theorem is far less 
restrictive than the condition for the factorizibility of the 
fermion determinant. Let us illustrate the general idea 
here by looking at the example of a two band spin 1/2 
model or a spin 3/2 model. In this case, we have fermion 
operators ipi p within one unit cell, where (3 = 1,2,3,4. 
Therefore, there are 16 fermion bilinears, of the form 

M 1 = 4,a M i^i^ where 1 = l.-,16. The 16 Af*j 
matrices can in general be expressed in a complete basis 
in terms of the product of s = 3/2 matrices Sf. 

I- 



S\ i = 1,2,3, 

£ijSiSj, a = 1, .., 5, = £ji, £jj = 0, 
t;fj k SiSjSk, L — 1, .., 7, £ij k — £j ik , £ iik = 0, 

(29) 

where £'s are fully symmetric, traceless tensors. If one in- 
sists on the factorizibility of the fermion determinant, one 
could only perform the HS decomposition in the density 
channel using the identity matrix /. However, since the 
£fjSiSj matrix contains an even power of spin matrices, 
it is also even under time reversal. HS decomposition in 
this channel does not lead to factorization of the fermion 
determinant, but according to our general theorem, it 
does lead to paired eigenvalues, and therefore, a positive 
fermion determinant. As we see from this non-trivial 
example, our method of T-invariant decomposition is in- 
deed more general and more powerful compared with the 
traditional method of factorization. We shall show the 
enlarged parameter space for QMC algorithm explicitly 
in the next section. 



IV. APPLICATION IN SPIN 3/2 AND n - 1/2 
HUBBARD MODEL 

In this section, we apply the method of T-invariant de- 
composition to the s = I model as an explicit example, 
and discuss the sign problem accordingly. After that, 
we generalize it to arbitrary fermionic Hubbard models 
with s = n — 1/2. These models are not of only aca- 
demic interests. In fact, the rapid progress in ultra-cold 
atomic systems provides an opportunity to study higher 
spin fermions. The simplest cases are the spin 3/2 atoms, 
such as 9 Be, 132 Cs, 135 Ba, 137 Ba atoms. Another impor- 
tant research direction is the trapped atoms in an optical 
lattice, formed by the standing wave laser beams, where 
the Hubbard model is a good approximation for these 
neutral atoms. 



A. The s=3/2 Hubbard model 

The spin 3/2 Hubbard model is defined as[| 

H = -t ^ { C ia C 3<r + h - c -} - (M + Mo) C ^ ClCT 
(ij),a i<? 

+Uo £ PfaWi) + U 2 p L(i)P2m(i), 

i i,m=±2,±l,0 

(30) 

with no — (Uq + 5C/2)/4. fi is fixed to be zero at half- 
filling on a bipartite lattice, to ensure the particle-hole 
(p-h) symmetry generated by the transformation Cj iCT — > 

(— ) l cJ CT . Because of the Pauli's exclusion principle, only 
on-site interactions in the total spin singlet (St = 0) and 
the quintet (St — 2) channels are allowed. PQ,P% m are 
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the singlet and quintet pairing operators defined by 
P o t (0(4o(i)) = ^={<| c l-|T<i c l-|}' 

P 2 (%) — 3C. i, -^lW = C ] 1 C ] -1 1 

b ) 2 ' 2 6 > 2 '» 2 

p 2 t ,_ 1 w = cU4-^ 4-2W = 4-i c l-i- ( 31 ) 



The s = 3/2 Hubbard model has an exact 50(5) or 
equivalently, Sp{A) symmetry, without any fine tuning 
of the parameters [J . This follows from the fact that 
singlet and quintet channel interactions can also be in- 
terpreted as 50(5) group's singlet and 5-vector represen- 
tations. When Uq — U2, the model has a larger symme- 
try, namely the S*t/(4) symmetry. The £*t/(4) symmetric 
Hubbard model has been extensively studied in the tran- 
sition metal oxides with double orbital degeneracy [l3| . 

To illustrate the T-invariant decomposition for this 
model, we first define the 4-component spinor 



ip(i) = (cz(i),ci(i),c _i(i),c _f (i))' 



(32) 



In this representation, we define five 4x4 Dirac r° (1 < 
a < 5) matrices to construct the Sp(4) or 5*0(5) algebra 
as 



r 1 = 



il 
-il 



• r 



2,3,4 



a 
-a 



,r 5 = 



-I 
-I 



where / and a are the 2x 2 unit and Pauli matri- 
ces. The ten S*0(5) generators are defined as r afc = 
-f[r a ,T b ](l <a,b< 5). Since the S*0(5) group is equiv- 
alent to the Sp(A) group, there exists a symplectic matrix 
i?, with the properties jlj 

R 2 = -1, i? f = R- 1 = l R=-R 

i?r a ir 1 = f r a , rt^r- 1 = r afc . (33) 

In our explicit representation, 



there are 16 bilinear fermionic operators, which can be 
classified into the scalar, vector, and anti-symmetric ten- 
sors (generators) of the 5*0(5) group as 



n(i) 
n a (i) 

L a b(i) 



1< a < 5 



1 



^(■Or&V7K«). l<a<^<5, (36) 



where n(i) is the particle number operator, n a (i) and 
L a b{i) represent the spin degrees of freedom. The ten 
S*0(5) generators are often conveniently denoted as 



/ Re7r x Re7r y Re7r z 











Uy 

S x 




V 



Im7r x 
Im7T y 
Ini7r z 




(37) 



Although the similar symbols are used in the S*0(5) al- 
gebra in the high T c cupratesQ, operators here have 
different physical meanings. 

These 16 fermion bi-linears are related through the 
Fierz identity 

E L ' h W+ E nl( l ) + j(n( l )-2) 2 = 5. (38) 



l<a<b<5 



Ka<5 



Defining the time reversal operator as T = RC, and using 
the properties of the R matrix given in Eq. I|33|) , it can 
be shown that n(i), n a {i) are even while L a b(i) is odd 
under T: 



TnT- 



Tn n T 



TL ab T 



J ab- 



(39) 



On the other hand, we can relate the above T-matrices 
with the usual spin SU{2) operators Ji, which form a 
subgroup of the S*0(5) group 



r = r x r 3 



-ia 2 
—i<j2 



(34) 



J-i- Jrf. it Z J lt 



Using the R matrix, the s = 3/2 Hubbard interaction 
can be written in an explicitly 5*0(5) symmetric fashion 
as 



tE^^ + tE^WW (35) 



where rf{i) = ip'(i)(R/2)ip^(i) is the singlet paring op- 
erator, and x a ^(i) = ift {i)(y a R/2)ip^ (i) are the polar 
forms of the quintet paring operators in Eq. (|3*T)) . 

In order to implement the method of T invariant 
decomposition, we first need to express the interac- 
tion terms in the particle-hole channel, rather than the 
particle-particle channel. In the particle-hole channel, 



Jz 



V3(-L 3 4 ± iL 2 i) 
—L23 + 2 L15 



(L12 ± iL 25 ) T i(Li3 ± ^35), 

(40) 



It is easy to check that the Si operators form the s = 3/2 
representation of the SU(2) algebra. While the above 
equation expresses the spin operators in terms of the T 
matrices of the S*0(5) algebra, the reverse can also be 
accomplished. The five T a matrices can actually be ex- 
pressed in terms of quadratic forms of the spin matrices: 



(41) 



where ff, is a rank-2 symmetric traceless tensor given 
Eq. (|2^|l . and discussed more explicitly in Ref. [l5| . 
The ten anti-symmetric tensor T ab matrices contain both 
the three linear (rank-1) Si and seven cubic symmetric 
traceless (rank-3) combination of SiSjSk operators, and 
they correspond to the second and the fourth rows of 
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Eq. J55Jl. Thus, the n(i) operator describes a particle- 
hole pair with total spin zero, the five n a (i) operators 
describe five particle-hole pair states with total spin two, 
and the ten L a t operators include the degenerate three 
spin-1 and seven spin-3 particle-hole pair states. From 
this point of view, the physical meaning of Eq. (|39|) and 
Eq. i|29[l becomes transparent: operators with even total 
spins are even under T, while operators with odd total 
spins are odd under T. 
Using the identities 



5 3 1 

(r a J R) Q/3 (i?r a ) 7d - = -SaySps - |r° 7 r^5 - ^r^r^, 



(42) 



and the Firez identity Eq. (|38[) . we can now express the 
s = 3/2 Hubbard model in the following form: 



Hj = 



i,l<a<5 



{|(n(z)-2) 2 + |n^)} (43) 



where 



9c = ~(3U + 5U 2 )/8 
g v = (U 2 -U )/2. 

B. Absence of the sign problem 



(44) 



After a series of transformations, we arrived at a form 
of the s = 3/2 Hubbard which is suitable for the T- 
invariant decomposition method. The interactions in Eq. 
(|43|l are fully expressed in T-invariant fermion operators 
in the particle-hole channel. When g c ,g v > 0, i.e, 



3/5U >U 2 >U , 



(45) 



the partition function can be expressed using the T- 
invariant decomposition as 

Z = J D^D^p exp{-J drii^ + H)^} 
= J DnJ Dn a exp[-^ ^ dr ][>(*, r) - 2) 2 

"Y f drJ2n 2 a ^r)}det{l + B}. (46) 

i,a 

. , „ , , . . o dT(Hff+Hl(T)) is obtained from 
the integration of fermion fields, n and n a are real HS 
bose fields. The time-dependent interaction -ffj(r) after 
the HS transformation is 

Hi{t) = -ff c ^ {T)ipi,<T (r)n(z, r) 

i 

~9v J2 4AT)r a a ^ a ,(T)n a (i, t). (47) 



3l t= " 5L i 


«2 









u„=u 2 



FIG. 2: Within the method of T-invariant decomposition, 
the shaded area marks the parameter region without the sign 
problem in the s — 3/2 Hubbard model at any doping level 
and lattice geometry. The fermion determinant can only be 
factorized along the ,5(7(4) line with Uo = U2 < 0, where 
the traditional algorithms can be applied without the sign 
problem. The sign problem along the line with Uo = Us > 
only disappears at the half-filling and on a bipartite lattice. 



We see that Hj(t) mixes the four spin components 
together, therefore the fermion determinant is factoriz- 
able if and only if g v = 0, which is the SU(A) line with 
Uo = U 2 < 0. We define the time reversal transformation 
T for the entire lattice as T = (J^ <&R(i))C. From Eq. 
(1391) we see that both terms are T-invariant: 



T{Hk + Hi)T~ 



H 



K 



Hr 



(48) 



and all other conditions of our theorem are met. There- 
fore, the minus sign is absent as long as Eq. (|45|l is 
satisfied. This is a much broader parameter range shown 
in Fig. |5J compared to the conventional factorizibility 
condition Uq — U 2 < 0. Our algorithm therefore enables 
us to study the s = 3/2 Hubbard model away from the 
SU(4) line. Our proof is valid for any filling level and lat- 
tice topology. In this parameter range, it is shown in Ref. 
1] that a number of interesting competing orders such as 
the staggered order of n a , singlet superconductivity, and 
the charge density wave can exist there. 

At half-filling and on a bipartite lattice where /i = 0, 
the sign problem also disappears along the SU (4) line at 
U = Uq = U 2 > 0. Similar to the spin 1/2 case, after 
performing a partial particle-hole transformation 



(49) 



while keeping c^a, Cji operators unchanged. The ki- 
netic energy part is invariant under the above transfor- 
mation, while the interaction part is change into Hi nt = 
2UJ2i -^15 W- It can b e decomposed using the imaginary 
number as 

Z = J T>Qexp{-2[/ J drQ 2 (i)|det|/ + s| 

where B = Te~ So dTH *+ H ^ with ff.( r ) = 
iUj^^Ur^^i^Q^T). Because T{iL 15 )T- 1 = 
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iLis, the det(/ + B) is positive definite. However, we 
did not succeed to generalize this at negative values of g c 
or g v away from the SU (4) line at half-filling. 

In practice, it is more efficient to sample with discrete 
HS transformation using two Ising-like fields 77, s for each 
quartic fermion term as in Rcf. 16] instead of using the 
continuous HS boson field. For any bilinear fermionic op- 
erator 0(i), the decomposition below has the numerical 
precision at the order of O(Ar) 4 as 



gArO(i,- 



t) 2 = ^2 lL e sr nV^90(i,T) + o(At 4 ), 



Ls=±l 



e -gArd(i,T) 2 _ lL e ™VlV&rgd(i,T) + ()(At 4 ), , 



l,8=±l 



where g > and 7; = 1 



2(3 -V6/)- Thc 



above proof for the positive definite of det(/ + S) applies 
equally well in this scheme. 

We only used the time reversal properties of the SO(5) 
algebra in the above proof; the exact SO(5) symmetry 
is useful for transforming the model expressed in the 
particle-particle channel to the particle-hole channel, but 
is not essential. A general anisotropic spin 3/2 lattice 
model is defined by 

H = - E {*ife+^to^+ |i - c '} 
+ E i h * n *M - " »(*)} + E { - f wo - 2 ) 2 

i.a i : a<b 

- f^(i) + ^L 2 o6 (i)}, (50) 

where t a is the spin dependent hopping amplitude, h a 
is the analogy of the Zeeman field coupling to the n a (i) 
field, g a and g a b are coupling constants in correspond- 
ing channels. When g c ,ga,gab are arbitrary positive 
interaction parameters, we can perform the same de- 
composition process as before. By using the fact that 
T(iL a f,)T _1 = iL a b, we again reach the positive defi- 
nite fermion determinant. This conclusion also holds for 
any valid representation of T matrices, with the rede- 
fined n(i),n a (i), L a b(i) and time reversal operations ac- 
cordingly. 



C. General higher spin Hubbard models 

We can generalize the results in the spin 3/2 case to 
the any fermionic system with spin s = n — | . The spin 
s = n — h Hubbard model can be written as 



H = -t^2 { C i* C 3° + h - c -} - (M + Mo) E 

{ij),cr ia 

+ E u j p jjS i ) p ^( i )> 



where J = 0,2,...,2n — 2 are the total spin of the particle- 
particle pairs, J z — 0, ±1, ...,±J. The pairing operators 



J,Jz 



are defined through the Clebsch-Gordan coefficient 



for two indistinguishable particles as 

a/3 



(52) 



The total spin of the particle-particle pair takes only 
even integer values so that Pauli principle is satisfied 
on every site. At half-filling and on a bipartite lat- 
tice, [i = ensures the particle-hole symmetry, and 
Mo = 1/ (2n) J2j(2J + i)Uj. 

The general strategy to implement the method of T- 
invariant decomposition is to first transform the inter- 
action terms originally expressed in the particle-particle 
channel to the particle-hole channel. In this case, we 
have fermion operators ipi p within one unit cell, where 
(3=1, .., (2s + 1). Therefore, there are (2s + l) 2 fermion 
bi-linears, of the form M 1 = ip\ a M^ 3 ijji t s, where I = 
1, (2s + l) 2 . The (2s + l) 2 M^ g matrices can in general 
be expressed in a complete basis in terms of the product 
of spin s matrices 



1, 

S\ » = 



1,2,3, 
a = 1, 



£ii,i 2 ,...ijSiiSi2 ' ' ' $i J) L — 1, •■, (4s + 1), 



(53) 



where £'s are fully symmetric, traceless tensors, satisfying 
£ L =£ L - (54) 

or any other permutation of indices, and 



Sii,ii, 







(55) 



Spherical harmonics can be used to explicitly construct 
these tensors [l7^. This decomposition is obviously com- 
plete, since 



(2s + l) 2 = l + 3 + 5 + ... + (4s + l) 



(56) 



According to the method of T-invariant decomposition, 
any negative interaction terms in the even spin channel 
like 1, £?jSiSj, ... or any positive interaction terms in the 
odd spin channel like Si, £,^ k SiSjSk, ■■■ can be simulated 
by our algorithm without the sign problem. 

In thc following, we shall illustrate this general proce- 
dure more explicitly for a special case of the higher spin 
Hubbard model where 



U 2 = U 4 



U2n-2 



u'. 



(57) 



(51) The generic higher spin Hubbard model only has the spin 
SU(2) symmetry for s 7^ |. However, under the above 
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condition, the higher spin Hubbard has the Sp(2n) sym- 
metry. When an additional condition, namely Uq = U' 
is imposed, the model has a larger, SU (2n) symmetry. 
In Appendix 1X1 an introduction to the Sp(2n) algebra is 
given. As shown there, the singlet pairing operator is also 
the singlet of the Sp(2n) group, while all other 2n 2 — n— 1 
pairing operators with J = 2, 4, 2n — 2 together form 
a representation for the Sp(2n) group. Thus we conclude 
that Eq. (|51|l is Sp(2n) symmetric if and only if coupling 
constants satisfy Eq. ijSTI) . For n=l and 2, the Sp{2n) 
symmetry is generic and does not need any fine-tuning. 
Actually, Sp(2) is isomorphic to SU(2), while 5p(4) is 
isomorphic to 50(5). This is consistent with our earlier 
finding that the s — 3/2 Hubbard model has the 50(5), 
or the 5p(4) symmetry without any conditions on the 
parameters 

To show the Sp(2n) symmetry explicitly, we can 
rewrite the Hamiltonian in Eq, (|51|l as 

H = -t^ {^ta^a + h.C.} -/J^4 ia ^,a 

(<i) i 
+ co ^(i^, a - nf + c 2 Y,^LY^,0) 2 , 

i i 

n + l TT 2n 2 -n-\ T . U -U' fro . 

The expression for Y a (l < a < 2n 2 — n — 1) is given in Ap- 



pendix^ where it is also shown that they are even under 
the time-reversal transformation. By the same reasoning 
before, we perform the HS decoupling in the above two 
channels. Then the sign problem is absent at both c\ and 
C2 are negative, i.e. 



Un < U' < - 



1 



2n 2 — n — 1 



U . 



(59) 



Because the Y a (1 < a < 2n 2 



- 1) are even under time- 



reversal transformation while the spin operators Si(i 
x, y, z) are odd, Y a can be expanded in the basis of Eq. 
B53J1 including all terms with even powers of spin matrix. 



V. BI-LAYER S=l/2 MODELS 

The spin 3/2 Hubbard model has a close relationship 
with a bi-layer model introduced by Scalapino, Zhang 
and Hanke(SZH)[T3. This model was constructed and 
extensively investigated because of the exact particle- 
particle channel 50(5) symmetry between the antifcr- 
romagnetism and superconductivi ty w hen the coupling 
constants satisfy a simple relation [is! flfll I2H. l2lj . The 
original SZH model was introduced on a two leg ladder, 
and it is straightforward to generalize it to a bi-layer sys- 
tem as 



H = -t|| { c L c 3<y + 4x d j> + h - c -} ~ t± { c la d ^ + h.c.) -fJ-Yl {™ C W + Ud ^} 

i 

+ V^(n c (i)-l)(n d (i)-l) + J^S itC -S itd , (60) 

i i 



where and eft are creation operators in the upper and 
lower layer respectively, £y and t± are the hopping am- 
plitude in the layer and cross the rung respectively, U 
is the onsite interaction, V and J are the charge and 
Heisenberg exchange interaction across the rung respec- 
tively. The SZH model is known to have an exact 50(5) 
symmetry when 

J = 4(U + V), fi = 0, (61) 

which unifies antiferromagnetism with superconductivity 
[Tl| . Remarkably, there exists another exact 50(5) sym- 
metry in the particle-hole channel when 

J = A{U-V), i± = 0, (62) 

and the symmetry is valid for all filling factors. We 
denote the former particle-particle 50(5) symmetry as 
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FIG. 3: The SZH model defined on a two-leg ladder segment 
of the double- layer spin 1/2 system. 



SO(5) pp and the later p-h 50(5) symmetry as SO(5) p h- 
The two 50(5) symmetric lines are shown in Fig. ^ In 
order to employ the method of T-invariant decomposi- 
tion, we adopt the view from the 50(5) symmetry in the 
particle-hole channel in this section. 

There are four single fermion states per unit cell 
in both the s — 3/2 Hubbard model and the s = 
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FIG. 4: Two 50(5) lines are shown in the SZH model as well 
as the QMC region without the sign problem for any filling 
(hatched area): g v \ > 0,g V 2 > and g c > 0. There is another 
region with V < (not shown). 

1/2 bi-layer model, a mapping between them can be 
established through ip^ = {c i ^,c i ^,c i _^,c i _^f r <-» 

(cj,T> c i,l' di,f) di,i) T ■ We denote the time reversal opera- 
tor defined in Eq. (|34|1 for the s = § system as T\, and 
the usual definition for s = 1/2 system as T%. T\ actu- 
ally is the combined operation of T2 and the interchange 
between the upper and lower layers: 

7i=(J,Qt 2 (63) 

The 16 p-h channel fermionic bilinear forms are 
mapped onto 



n(i) 


= 4<T C i°- + 4a d i° 




= -i{d\ a c ia -h.c.)/2 


1*5(1) 


= (4 a c ia + h.c)/2 


"2,3,4(0 




Re7r(i) 


= 4 a (^)^d i0 + h.c. 


Im7?(i) 




Si 


= c L(|)«/3Ci/9 + d\ a {~) a pd ip 


Q 


= (n c (i) - n d (i))/2, 



(64) 



where 71^5 are the singlet rung current and rung bond or- 
der parameters respectively, and Im7? and Re7? are their 
triplet counterpart; 71-2.3,4 are the rung Neel order pa- 
rameter, and S is the total rung spin, Q is charge density 
wave order parameter, n, are even under the T\ 

transformation and the others are odd. In contrast, n, 
n.5, Imjf and Q are even under the usual definition T2 
and the others are odd. 



The general SZH model can be mapped into an 
anisotropic SO(5) model in the following form: 

(ij) ' » 

+ E { - 1 wo - 2 ) 2 - ir(»?(o + »2(0) - ^ 

i 

x (n|(i)+n2(0+n|(0), (65) 
with 

Ag c = \j-U- 3V, 4g vl = ^J-U + V, 
4^2 = 'l + U-V. (66) 

The particle-hole channel SO(5) p h symmetry is restored 
at g v i = g V 2 and t j_ = 0, i.e., when the conditions of 
Eq. I|62|) are satisfied. At this point, the SZH model 
expressed in Eq. i|65|) takes exactly the same form as 
the s = 3/2 Hubbard model expressed in Eq. H43fl . The 
equivalence between the two models is therefore rigor- 
ously established. 

For the general SZH model, the interactions can be 
expressed purely in terms of the fermion bi-linears which 
are invariant under T\ transformation from Eq. 16511 . 
by virtue of Eq. I|39|) . We perform the T\ invari- 
ant decomposition of the interactions in the region of 
9c,9vi,9v2 > 0, i.e. 

(y+v>u>-y+v 

\\j>u + w ' (b7j 

as shown in Fig. 0] then the sign problem is absent. In 
this region with positive g v i^ and g c , we expect the com- 
peting orders of the five- vector channel and the supercon- 
ductivity, i.e. the antifcrromagnctism, staggered current 
and the rung-singlet superconductivity, which can be in- 
vestigated systematically with high numerical accuracy, 
at and away from the half-filling. 

The above algorithm has been applied to demonstrate 
the existence of 2-dimensional staggered current phase 
conclusively at half-filling with t\\ = l,t± = 0.1, U = 
0, V = 0.5, J = 2 [iJJ]. The current pattern is illustrated 
in Fig. with staggered inter-layer currents (SIC) be- 
tween the bi-layers and alternating source to drain cur- 
rents within the bi-layers. Viewed from the top, this 
current pattern has a s-wave symmetry. While the D- 
density wave |22| currents are divergence free within the 
layer, the SIC current is curl free within the layer. These 
two patterns can be considered as dual to each other in 
two dimensions. As far as in our knowledge, this is the 
first time a current carrying ground state has been con- 
clusively demonstrated in a 2D system. 

The mapping between the SZH model and the s = 3/2 
model are not unique. More generally the SZH model 
can be written as 

(ij) i 1 
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FIG. 5: (a) Sketch of a staggered interlayer current phase 



a]_Sl 
10]. 



from Ref. |lf| . For clarity, we do not show the bottom layer 
current. (b)Top view of the bi-layer current, (c) Sketch of the 
D-density wave current pattern for comparison. 



+ E{- - 2 ) 2 - T (n?(<) + ^(0) - 9 f 

i 

x («*(») + nl(i) + n\{i)) + 9 fQ 2 {i) + 9 -f{S{i) ■ S{i)) 
+ ^(Re7f(i)-Re7r(i)+Im7r(i)-Im7f(i)). (68) 

Only three out of the six coupling constants are indepen- 
dent, as shown here in the correspondence to the U, V, J 
parameters 



U = -4g c + 3.^2 + g t i - 3gt2 
V = -4g c + g v i - 9tx - 3g f 3 
J = A(g v i + g v2 + g t 2 + gti)- 



(69) 



If the gt\i gt2,9t3 are set to zero, it returns to Eq. (|65ll . 
For any given values for U, V, J, if we can find a set of 
value of g c , g v i, g V 2, gti, gt2, gt3 > 0, then we can per- 
form the HS transformation keeping the invariance un- 
der the T\ operation and arrive at the absence of the 
sign problem regardless the doping and lattice topology. 
This general decoupling scheme extends the valid param- 
eter region in Eq. I|67l) . On the other hand, we can also 
consider to perform HS decoupling with the invariance 
under the usual definition of the time reversal operation 
T2. After setting g v \, g t 3 = 0, we have the condition that 
9c > 0,g V 2 < 0,541 < 0, gt2 > 0. This decoupling scheme 
based T\ also enlarges the region of Eq. lJS7)l. For ex- 
ample, the usual bilayer negative U Hubbard model with 
U < 0,V — J — is out of that region. Nevertheless, we 
can still show the absence of the sign problem by setting 
g c = — gn /A = — U/2 > and all other parameters zero. 

In contrast, the conventional algorithm based on the 
factorization of the fermion determinant works only at 



either g v \ — g V 2 — 0,g c < or the usual negative U 
Hubbard model U < 0, V = J = 0. This parameter 
set is included in the above HS decomposition schemes 
respecting cither the T2 or the T\ time reversal symmetry. 
We therefore see the significant improvement provided by 
the method of the T-invariant decomposition. 
VI. MODELS WITH BOND INTERACTIONS 

So far, the models we considered only have on-site in- 
teractions. In this section, we will generalize them to 
include interactions defined on the bond. Such models 
can have many exotic phases. 

We first consider the following general single layer spin 
1/2 Hamiltonian with bond interactions 

H = -t E (4* c 3° + h.c.) - n^n{i) 

(ij),* i 

{ _ i^ Mi jMi , + ^pv\^ 



E 

(ij) 
gtbd 



Y 



M ij = 4*(|W c j73 + h.c, Nij = i{4j^) a pCj P - h.c.}, 



My = c\ a Cj a + h.c, = i{cl a Cj a - h.c.}, 



(70) 



where My and iVy are the singlet bond and current op- 
erators on the bond (ij), Mij and iVy are their triplet 
counterparts, g s bd, g SC ur, gtbd and g tC ur are the coupling 
constants in the corresponding channels. The sites i and 
j forming the bond (ij) are not necessary nearest neigh- 
bors, but can be at arbitrary distance apart. Under the 
time reversal transformation T, My and Nij are even 
while Mij and JVy are odd. 

These four interactions are not independent, and can 
be reorganized into 

H int = E { - J c( c lT c !i c il c iT + h - c -) 

(ij) 

+ V(n(i) - l)(n(j) - 1) + J s S(t) ■ S(j), } 

Jc = ^{fjsbd + 9scur) + %{fjtbd + gtcur)-, 
Qsbd Qscur 3 



v = 



-^{gtbd — gtcur), 



Js — 2(g S bd — gscur) + (gtbd — gtcur)- 



(71) 



The J c term is the pair hopping, V is the charge inter- 
action between site i and j, and J s is the Heisenberg 
exchange. When all of g sbd , g scur , g t bd, gtcur are positive, 
we perform the HS decomposition in each channel respec- 
tively as 
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Z = DMDMDNDN cxp { - / dr ff.wAfg (r) + gaBur N§ (r) + g t bdM% (r) + gtcurN* (r) 2 } 

x det{I + B}, (72) 
where I + B = I + Te Jo dTHK+Hl ( T \ Hj(t) after the HS decoupling is given by 

Hi{t) = - g s bdMjj (t) (c| -ct Cj> + /i.e.) +«^Sscuri%(T)i(ct iCT Cj )(T - /i.e.) 

(y) ~ <»i) 





C 



D 



FIG. 6: Four possible density-wave phases can be simu- 
lated without sign problem. A) singlet spin-Peierls (p-density 
wave), B) Triplet d x 2_ y 2 density wave, C) singlet d xy density- 
wave, D) Triplet diagonal current. 



Therefore, Hi and I + B are even under the time- reversal 
transformation and the sign problem is absent. 

The valid parameter region for the above algorithm is 
very general as long as all g sbd , g scur , g tbd , g tcur > 0. As a 
result, V and J can be either positive or negative while 
J c has to be positive. Many interesting competing orders 
are supported in this parameter region. For example, 
various density-wave states exist on a square lattice near 
half-filling 23| as shown in Fig. HJ1 With g s bd>9tcur > 0, 
the above algorithm provide a good opportunity to study 
the singlet bond and the triplet current order parame- 
ters formed by My and 2Vy , while its not good for study 
the singlet current and the triplet bond order param- 
eters formed by Nij and My- because gtbd, gscur > 0. 
After setting g t bd = gscur — 0, for the bond (ij) connect- 
ing the nearest sites, the g sb d term favors the p-density 
wave (spin-Peierls) phase, and the gtcur term favors the 
triplet channel d a; 2_ v 2-density wave. The latter order is 



recently proposed as the origin as the pseudogap in the 
high T c cuprates|2^|. For the bond interaction between 
the next nearest bond, i.e. the diagonal bond, the g s bd 
term leads to the singlet d xy order, and gtcur term leads 
to the triplet diagonal current order. The triplet diagonal 
current phase was studied in the two-leg ladder system 
using the bosonization method in Ref. |25j and also un- 
der the name of the triplet F-density wave in Ref. j2|J. 

When the Fermi surface nesting effect is not impor- 
tant either at large doping or in the non-bipartite lat- 
tice, the gtcur term can lead to the F" channel of the 
Landau-Pomeranchuk instability on the Fermi surface, 
which was studied recently in the continuum model in 
Ref. |27|. After the symmetry breaking, two possible 
phases are named as a and f3 phases in analogy to the 
A and B phases in the triplet p-wave channel superfluid 
phase in 3 He as shown in Fig. [7| The a-phase was stud- 
ied by Hirsch 28, 29j under the name of spin-split phase 
on the lattice system with an opposite anisotropic Fermi 
surface distortions for two spin components. In contrast, 
the fermi surface distortion is isotropic and a spin-orbit 
coupling is dynamically generated in the (3 phase. The 
two single particle bands arc characterized by the hclici- 
ties. It would be interesting to study these exotic phases 
in our version of the QMC algorithm free of the sign 
problems. 

Bond interactions can also be added into the spin 3/2 
Hubbard model of Eq. (|30f) as 



Hi 



bond 



(ij) 



i^ Mij M ij + ^N ijNij 



9vbd , , a nja gvs 



Sp_yvbd M a M 



2 —tj"^3 ' 2 ~ Ni ° Nlj 

9tbd.rabj.rab 9tcur 



a<b 



2 13 



$N$}, (74) 



Mij = ipjipj + h.c, = i{iplipj - h.c.}, 



M a 



M; 



ah 



+ h.c, N« = i{4—^ - h.c.}, 
+ h-c, N$ = i{4^-1>, ~ h.c.}, 
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a-phase fj-phase 

FIG. 7: The fermi surface instability in the F£ channel, 
with dashed lines marking the fermi surface before symme- 
try breaking. In the a-phase, the anisotropic fermi surface 
distortion appears for two spin components. In the /3-phase, 
spin-orbital coupling is generated dynamically and two fermi 
surfaces are characterized by helicity. 



be simulated, and some novel and exotic ground states 
have been firmly established. 

We conclude this paper with an optimistic outlook. 
Even though our method can only be applied presently 
to models with definite constraints among the interaction 
parameters, we believe that the deep symmetry connec- 
tions revealed in this work could guide us in future works, 
and might eventually lead to the complete elimination of 
the sign problem. 

Note added in proof: 
After the paper has been accepted, we learned that a 
similar version of the theorem of T-invariant decomposi- 
tion had been discussed in the context of lattice gauge 
theory [3(j • However, they did not consider the case that 
I+B is not diagonalizable. Our proof is valid regardless 
of whether I+B is diagonalizable or not, thus is more 
complete." 



(75) 

where My and Nij are the singlet bond and current op- 
erators on the bond (ij), M§,N§, M?j>, N$ are their 
5-vector and 10-tensor channel counterparts respectively, 
and g sbd , g vbd , gtbd, 9scur, gvcur, gtcur are the coupling con- 
stants in the corresponding channels. Again the site 
i and j forming the bond (ij) can be at an arbitrary 
distance apart. The bond interactions can be decou- 
pled by introducing the HS field in each channel respec- 
tively. Following the same reasoning as the case of the 
spin 1/2, the bond-interactions keep the fermion deter- 
minant positive definite provided all these coupling con- 
stant non-negative. Similarly, with g s bd, gvbd, 9tcur > 0, 
the algorithm can be applied to study the singlet, quin- 
tet bond orders and the 10-fold current order, while with 
g S cur,gvcur,gtbd > 0, it is not useful for study applied to 
study the singlet, quintet current orders and the 10-fold 
bond order. 



VII. CONCLUSION 

The sign problem of the fermionic QMC algorithm 
is one of the most important problems in theoretical 
physics. Its solution would practically give an universal 
computational method to solve models with strong cor- 
relations. The rigorous theorem established in this work 
shows that the minus sign problem can be eliminated for 
a much wider class of models than before, in which the 
fermion matrix is invariant under an anti-unitary sym- 
metry similar to the time reversal symmetry in quantum 
mechanics. The method of T-invariant decomposition 
does not only provide a deep connection between the sign 
problem and the time reversal symmetry, it also leads to 
practical algorithms which can be applied to many inter- 
esting models with strong correlations. Using this algo- 
rithm, a new class of models with strong correlation can 
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APPENDIX A: Sp(2n) ALGEBRA IN THE SPIN 
s = n - 1/2 FERMION SYSTEM 

We give a brief introduction to the Sp(2n) algebra here. 
The 2n-dimensional Hilbert space on each site can be ar- 
ranged as a direct product between a n-dimcnsional and 
a 2-dimensional space. The complete basis of eigenstates 
of S z are labeled in the sequence of |1) = |n — |), |2) = 

I -" + §>. \n) = l^r 11 ), and |1> = I - n + 3>» I 2 ) = 
|n — |), |n) = l^-jf— ). The Sp(2n) spinor is defined as 

V> = (c„_i,c_ n+ i,c_„ + 3,c„_3,....) T . (Al) 

Group elements of Sp(2n) include any 2n x 2n uni- 
tary matrix U satisfying U T RU — R or equivalently 
R- 1 UR=U* [3ll with the i?-matrix 

R = I n ®(-io- 2 )- (A2) 

The R matrix is a straightforward generalization of the 
R = —io-2 in the spin 1/2 case, which also satisfies 
R T = R- 1 = i? f = -R. Clearly, the Sp(2n) group is 
a subgroup of the SU(2n) group defined in the 2n di- 
mensional space. 

In the particle-hole channel, there are 4n 2 independent 
bilinear operators as V'aV' p{ a — 1, ...,2n,/3 — 1, ...,2n). 
Among them, the particle density operator n = ifi^/ipa is 
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a singlet under both the SU (2n) and the Sp(2n) group. 
The time reversal transformation T = CR is defined as 
usual, and it satisfies T 2 — — 1. The other An 2 — 1 bilinear 
operators form the generators (adjoint representation) for 
the SU(2n) group. They can be decomposed into two 
classes according to their transformation properties un- 
der the T operation. The first class contains n(2n + 1) 
elements which forms the generators of the Sp{2n) group 
as denoted as ip^X^gippib = 1 ~ 2n 2 + n). X b can be 
expressed in terms a direct product between the SU (n) 
and SU{2) generators. We define the SU(n) generators 
as 

{M^)ik = ^(SuSjk + StkSji) (l<i<j< n), 

(Mg^ik = ^-(SaSjk - SikSji) {l<i<j< n), 

w(3) diag(l,...,l,-(j-l),0,..,0) 

Mi = (2 < 7 < n), 

(A3) 

where M§ } , Mfp , Mj 3) arc the n x n dimensional gener- 
alization of the SU (2) Pauli matrices cr x ,y,z respectively. 
Counting the numbers of SU{2n) generators, there are 
n(n — l)/2 real symmetric M^'s, n(n — l)/2 imaginary 
anti-symmetric Mfj's, and n — 1 real diagonal M^-'s. 
Then the Sp(2n) generators X h can be expressed as 

Mlf®I 2 , M^Qff, Mf ] ®a, I n ® a, (A4) 

Because Af^'s and M, (3), s are real and M^ 2 ' 's are purely 
imaginary, the Sp(2n) generators are odd under the time 
reversal: T~ 1 X b T = —X b . The second class bilinears 



V'a^aflV'jS nav e 2n 2 — n — 1 elements. y a (a = 1, 2n 2 — 
n — 1) are given by 

M^Wj Afj^g/a, A/j 3) <g>/ 2 , (A5) 

which are even under the time reversal: T _1 Y a T = Y a . 

These An 2 bilinear operators are not independent of 
each other, but are related by the Fierz identity. The 
total Hilbert space for one site has the dimension of 2 2n , 
which can be decomposed into subspaces with different 
particle number r(0 < r < 2n). Each of them form 
the totally anti-symmetric representation of the SU(2n) 
group V. The Casimir value of the SU(2n) group in 
such representations are r(2n + l)(2n — r)/(2n). Thus 
we arrive the Fierz identity for the spin n — \ system as 

b a 

+ ^^(^^«-«) 2 = g ' (A6) 

The onsite pairing operators can be easily formed 
by using the R matrix. Due to the Pauli's exclusion 
principle, the total spin for a s-wave pair can only be 
0, 2, 2n— 2. The singlet pair operator is also the Sp(2n) 
singlet operator. It can be written as tp^Rapipt, which 
was studied extensively in a Sp( 2n) generalization of the 
Heisenberg antiferromagnet [32( . The other 2n 2 — n — 1 
pairing operators with total spin 2, 4, 2n — 2 together 
form a representation of Sp(2n) as ^(^(RY^upTlrJi). 
When all the interaction parameters are equal, these 
n(2n— 1) pairing operators together form anti-symmetric 
representation of SU{2n) of V (r=2). 
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